# Figure_2.R
# (formerly normal_normal_divergence_code.R)
#
# 2009 July 21
# John G. Bullock
#
# Creates PDF of Figure 2 in "Partisan Bias and the Bayesian Ideal
# in the Study of Public Opinion," Journal of Politics 71 (July): 
# 1109-24, http://dx.doi.org/10.1017/S0022381609090914.
#
# sessionInfo() output is at the bottom of this file.

if (any(regexpr('Hmisc', search())>-1)) { detach(package:Hmisc) }
library(grid)
library(lattice)
set.seed(1977) # to get different simulation results, change the seed or delete this line


# GENERATE THE DATA
nrows <- 5
ncols <- 4
panel.periods <- rep( 6, nrows*ncols) # for how many stages will we see updating in each panel?
simsperpanel  <- 030                  # number of simulations plotted in each panel

D.prior.means     <- rep(2, nrows*ncols)
D.prior.vars      <- rep(c(.5, 1, 1.5, 2), nrows)
R.prior.means     <- rep(1, nrows*ncols)
R.prior.vars      <- rep(1, nrows*ncols)
mu      <- 2
x.vars  <- rep(c(1, 5, 10, 15, 20), 1, each=ncols)
x.means <- D.post.means <- R.post.means <- vector("list", nrows*ncols) # set up an empty list with one element for each panel

Bayes.update <- function(prior.mean, prior.variance, m, m.variance) {
  post.means <- rep(NA,length(m))
  post.variances <- rep(NA, length(m))
  
  post.means[1]     <- prior.mean*((1/prior.variance)/((1/prior.variance)+(1/m.variance))) + m[1]*((1/m.variance)/((1/prior.variance)+(1/m.variance)))
  post.variances[1] <- 1/(1/prior.variance + 1/m.variance)
  if (length(post.means)>1) {
        for (i in 2:length(post.means)) {
              post.means[i]     <- post.means[i-1]*((1/post.variances[i-1])/((1/post.variances[i-1])+(1/m.variance))) + m[i]*((1/m.variance)/((1/post.variances[i-1])+(1/m.variance)))
              post.variances[i] <- 1/(1/post.variances[i-1] + 1/m.variance)
        }
  }
  post <- cbind(post.means, post.variances)
  return(post)
}

for (i in 1:length(x.means)) {
  x.means[[i]] <- D.post.means[[i]] <- R.post.means[[i]] <- vector("list", simsperpanel)
  for (j in 1:length(x.means[[i]])) {
    x.means[[i]][[j]]      <- rnorm(n=panel.periods[i], mean=mu, sd=sqrt(x.vars[i]))
    D.post.means[[i]][[j]] <- c(D.prior.means[i], Bayes.update(prior.mean=D.prior.means[i], prior.variance=D.prior.vars[i], m=x.means[[i]][[j]], m.variance=x.vars[i])[, "post.means"])
    R.post.means[[i]][[j]] <- c(R.prior.means[i], Bayes.update(prior.mean=R.prior.means[i], prior.variance=R.prior.vars[i], m=x.means[[i]][[j]], m.variance=x.vars[i])[, "post.means"])    
  }
}


# CREATE THE DATA FRAME
normal.normal.divergence.data.fixed <- expand.grid(stage=0:panel.periods[1], simulation=1:simsperpanel, panel=1:20) # stage 0 gives the priors
normal.normal.divergence.data.fixed$mu.D <- unlist(D.post.means)
normal.normal.divergence.data.fixed$mu.R <- unlist(R.post.means)
normal.normal.divergence.data.fixed$diff <- abs(unlist(D.post.means)-unlist(R.post.means))
normal.normal.divergence.data.fixed$diff.from.prev <- c(NA, diff(normal.normal.divergence.data.fixed$diff))
  normal.normal.divergence.data.fixed$diff.from.prev[normal.normal.divergence.data.fixed$stage==0] <- NA
normal.normal.divergence.data.fixed$var.D.prior <- rep(rep(D.prior.vars, each=panel.periods[1]+1), each=simsperpanel)  
normal.normal.divergence.data.fixed$var.R.prior <- rep(rep(R.prior.vars, each=panel.periods[1]+1), each=simsperpanel)  
normal.normal.divergence.data.fixed$x.vars      <- rep(rep(x.vars,       each=panel.periods[1]+1), each=simsperpanel)  


# FIND MINIMUM AND MAXIMUM RANGE OF VALUES IN EACH PANEL -- FOR INFORMING READERS ABOUT THE RANGE ON THE Y-AXIS
# tmpmin <- rbind(  tapply(normal.normal.divergence.data.fixed$mu.D, normal.normal.divergence.data.fixed$panel, min)
#                 , tapply(normal.normal.divergence.data.fixed$mu.R, normal.normal.divergence.data.fixed$panel, min))
# tmpmin <- apply(tmpmin, 2, min)
# tmpmax <- rbind(  tapply(normal.normal.divergence.data.fixed$mu.D, normal.normal.divergence.data.fixed$panel, max)
#                 , tapply(normal.normal.divergence.data.fixed$mu.R, normal.normal.divergence.data.fixed$panel, max))
# tmpmax <- apply(tmpmax, 2, max)
# which((tmpmax-tmpmin)==min(tmpmax-tmpmin))
# which((tmpmax-tmpmin)==max(tmpmax-tmpmin))
# rbind(tmpmin,tmpmax) # range from min(min(D), min(R)) to max(max(D), max(R))
# diffmin <- tapply(normal.normal.divergence.data.fixed$diff, normal.normal.divergence.data.fixed$panel, min)
# diffmax <- tapply(normal.normal.divergence.data.fixed$diff, normal.normal.divergence.data.fixed$panel, max)
# rbind(diffmin, diffmax)


# WITHIN-STAGE DIFFERENCES BETWEEN UPDATERS
# tapply(normal.normal.divergence.data.fixed$diff, normal.normal.divergence.data.fixed$panel, max)   # largest within-stage differences
# normal.normal.divergence.data.fixed$diff[normal.normal.divergence.data.fixed$stage==panel.periods] # differences at final stage
# tapply(  normal.normal.divergence.data.fixed$diff.from.prev
#        , normal.normal.divergence.data.fixed$panel, function (x) any(x>0))                         # which panels show stage-to-stage divergence?



# PRELIMINARIES FOR PRINTING THE FIGURE
dir.output      <- '.' # where the figures go
filename.output <- 'Figure_2.pdf'
PDFtitle        <- 'Divergence and convergence under Bayesian learning about an unchanging condition'
PSwidth         <- 9.5
PSheight        <- 12.0
postscript.bg   <- 'transparent'
panelwidth      <- list(1.090, "inches")       # these dimensions allow 4 columns on 8.5" paper
panelheight     <- list(0.86875, "inches") 
plotlinealpha   <- .725                        # transparency (1 = opaque)
plotlinecolor   <- rgb(0, 0, 0, plotlinealpha) # first three params are gray (lower=darker)
plotlinewidth   <- 1
panelbordercol  <- 'black'
panelnumbercol  <- grey(.2)                    # color of number of each panel; lower values are darker
panelborderwidth <- .3
panelLayout     <- c(ncols, nrows)             # columns, then rows
xbetween        <- 1.105                       # space between columns
ybetween        <- 1.3175                      # space between rows
panel.xlim      <- NULL                        
panel.ylim      <- c(0, 1.75)  
striptextsize   <-  .75                        # cex
xaxtextsize     <-  .8*.87                     # cex
yaxtextsize     <-  .8*.87                     # cex
axisticksize    <- .4
xaxs            <- list(draw=T, labels=0:panel.periods[1], at=0:panel.periods[1], tck=c(axisticksize, 0), col=panelbordercol, cex=xaxtextsize,
                        alternating=1, relation='free', axs='i') # axs='i' means that there is no padding around the xlimits
yaxs            <- list(draw=T
                        , labels=c(0, '.5', '1', '1.5'), at=c(0, .5, 1, 1.5)
                        , tck=c(axisticksize, 0)
                        , col=panelbordercol
                        , cex=yaxtextsize
                        , alternating=1
                        , relation='free' # but the ylim will constrain all panels to same height
                        , rot=90 #
                        , axs='i'
                        )


# FUNCTION TO STRIP BAD CHARACTERS FROM BEGINNING OF STRING
stripinit <- function(x,bad=" ") { # strips bad from beginning of string
  bad <- paste(bad,"*",sep="")
  sub(bad,'',x)
}


# PANEL FUNCTION
normal.normal.divergence.panel <- function(...) { 
  trellis.par.set("clip", list(panel="off", strip="off")) # for permitting panel.axis(outside=T, ...)
  panel.xyplot(...)

  # y-axis
  trellis.par.set("clip", list(panel="off", strip="off")) # for permitting panel.axis(outside=T, ...)  
  grid.text('difference of means', x=-.240, y=.48, hjust=.5, vjust=.5, gp=gpar(font=1, cex=.705, col=panelbordercol), rot=90)    
 
  # per-panel x-axis labels
  grid.text('time', x=.50, y=-.33, just='center', gp=gpar(font=1, cex=.72, col=panelbordercol))

  # labels above each column and beside each row
  # note: expression(rho) produces \sigma
  var.D.prior <- stripinit(normal.normal.divergence.data.fixed$var.D.prior[normal.normal.divergence.data.fixed$panel==panel.number()][1], 0)
  var.R.prior <- stripinit(normal.normal.divergence.data.fixed$var.R.prior[normal.normal.divergence.data.fixed$panel==panel.number()][1], 0)  
  x.var       <- stripinit(normal.normal.divergence.data.fixed$x.vars[normal.normal.divergence.data.fixed$panel==panel.number()][1], 0)

  if (panel.number()%in%1:4) {
    grid.text(substitute(sigma[D0]^2==var.D.prior), x=.5, y=1.25, just='center', gp=gpar(font=1, cex=.8, col=panelbordercol))
    if (panel.number()==1) {
      grid.lines(x=c( .005, 5.305), y=c(1.10,  1.100), gp=gpar(lwd=0.75, col=panelbordercol, lineend='square')) # horizontal line
      grid.lines(x=c(-.400, -.400), y=c( .98, -6.455), gp=gpar(lwd=0.75, col=panelbordercol, lineend='square')) # vertical line
    }    
  }

  if (panel.number()%in%c(1,5,9,13,17)) { grid.text(substitute(sigma[x]^2==x.var), x=-.50, y=.5,  vjust=.5, hjust=.5, rot=90, gp=gpar(font=1, cex=.8, col=panelbordercol)) }

  # miscellaneous labels
  grid.text(panel.number(), x=.965, y=.875, just='right', gp=gpar(font=1, cex=.80, col=panelnumbercol)) # lower numbers are darker for grey()
  if (panel.number()==18) { grid.text(paste(Sys.time(), '; mu = ', mu, '; D.prior = ', D.prior.means[1], '; R.prior = ', R.prior.means[1], '; simsperpanel = ', simsperpanel, '; alpha = ', plotlinealpha, '; lwd = ', plotlinewidth, sep=''), y=-1, just='center', gp=gpar(font=2)) }  
}


# MAKE THE PDF FILE
setwd(dir.output)
pdf(file=filename.output, width=PSwidth, height=PSheight, paper="special" 
           , title=PDFtitle
           , bg=postscript.bg)
trellis.par.set("dot.symbol", list(alpha=1, cex=0, col="black", font=1, pch=16))                      # black dot (not grey) to mark the mean / cex=0 for no dot; .5 is normal
trellis.par.set("axis.line", list(alpha=1, col=panelbordercol, lty=1, lwd=panelborderwidth))          # to eliminate lattice panel border, set col="white"
trellis.par.set(superpose.line=list(col=plotlinecolor, lty=c(1,1), lwd=plotlinewidth))      
normal.normal.divergence.plot <- xyplot(diff~stage | panel, data=normal.normal.divergence.data.fixed, type="l", 
                   , groups=simulation 
                   , xlab='', ylab=''                    
                   , ylim=panel.ylim
                   , main=''
                   , panel=normal.normal.divergence.panel
                   , layout=panelLayout               # x columns, y rows
                   , as.table=TRUE                    # plot panels left to right, top to bottom
                   , between=list(y=ybetween, x=xbetween)                   
                   , strip=F
                   , scales=list(x=xaxs, y=yaxs)
                   )
print(normal.normal.divergence.plot, panel.width=panelwidth, panel.height=panelheight) # more=TRUE
dev.off()
shell.exec(paste(getwd(), '/', filename.output, sep=''))   

# > session.info
# R version 2.6.1 (2007-11-26) 
# i386-pc-mingw32 
#
# locale:
# 		 LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
#
# attached base packages:
# 	 	 [1] grid      stats     graphics  grDevices utils     datasets  tcltk     methods   base     
#
# other attached packages:
#		 [1] lattice_0.17-2 car_1.1-2    
